#ifndef _NFEM_SF_H
#define _NFEM_SF_H

#include "exception"
#include "stdexcept"

#ifndef NFEM_SF_ZERO
#define NFEM_SF_ZERO 1e-10
#endif // NFEM_SF_ZERO

/**
 * @brief      Calculate the Jacoby value of a 2x2 matrix
 *
 * @param[in]  x00   Value of element (0, 0)
 * @param[in]  x01   Value of element (0, 1)
 * @param[in]  x10   Value of element (1, 0)
 * @param[in]  x11   Value of element (1, 1)
 *
 * @return     The Jacoby value
 */
double jacoby2(double x00, double x01, double x10, double x11);

/**
 * @brief      Calculate the Jacoby value of a 3x3 matrix
 *
 * @param[in]  x00   Value of element (0, 0)
 * @param[in]  x01   Value of element (0, 1)
 * @param[in]  x02   Value of element (0, 2)
 * @param[in]  x10   Value of element (1, 0)
 * @param[in]  x11   Value of element (1, 1)
 * @param[in]  x12   Value of element (1, 2)
 * @param[in]  x20   Value of element (2, 0)
 * @param[in]  x21   Value of element (2, 1)
 * @param[in]  x22   Value of element (2, 2)
 *
 * @return     The Jacoby value
 */
double jacoby3(double x00, double x01, double x02, double x10, double x11, double x12, 
	double x20, double x21, double x22);

// Space of node based finite element method shape (basic) function
namespace nfem_sf
{

/**
 * @brief     Enum type of calculation target
 */
enum valtype_e
{
	Value,
	Dx,
	Dy,
	Dz,
};

/**
 * @brief      The 1st order (linear) shape functions
 */
class linear
{
public:
	linear(){}
	virtual ~linear(){}

	/**
	 * @brief      Check to see if the input parameters of a 1-D line element is valid.
	 *
	 * @param[in]  x     independent variable
	 * @param[in]  x1    low bound of the shape function
	 * @param[in]  x2    high bound of the shape function
	 *
	 * @return     true if the inputs are valid.
	 */
	bool line_element_is_valid(double x, double x1, double x2);

	/**
	 * @brief      Shape function defined on a 1D line.
	 * 
	 * This function has two nodes on the left and right ends, respectively.
	 * index for the left node function is 0. And the node's index is 1.
	 * 0-----------1-----> x
	 *
	 * @param[in]  x     independent variable
	 * @param[in]  x1    low bound of the shape function
	 * @param[in]  x2    high bound of the shape function
	 * @param[in]  idx   Node index
	 * @param[in]  vt    Calculation target. Could be Value or Dx
	 *
	 * @return     Value of the shape function
	 */
	double line(double x, double x1, double x2, unsigned int idx, valtype_e vt);

	/**
	 * @brief      Check to see if the input parameters of a 2-D quadratic element is valid.
	 *
	 * @param[in]  x     independent variable x
	 * @param[in]  y     independent variable y
	 * @param[in]  x1    low bound of x
	 * @param[in]  x2    high bound of x
	 * @param[in]  y1    low bound of y
	 * @param[in]  y2    high bound of y
	 *
	 * @return     true if the inputs are valid.
	 */
	bool quad_element_is_valid(double x, double y, double x1, double x2, double y1, double y2);

	/**
	 * @brief      Shape function defined on a 2D rectangle.
	 * 
	 * This function has four nodes on the corners, respectively.
	 * index of the nodes are as show
	 * 
	 * y
	 * |
	 * |
	 * 3----------2
	 * |          |
	 * |          |
	 * |          |
	 * |          |
	 * 0----------1--->x
	 *
	 * @param[in]  x     independent variable x
	 * @param[in]  y     independent variable y
	 * @param[in]  x1    low bound of x
	 * @param[in]  x2    high bound of x
	 * @param[in]  y1    low bound of y
	 * @param[in]  y2    high bound of y
	 * @param[in]  idx   Node index
	 * @param[in]  vt    Calculation target. Could be Value, Dx or Dy
	 *
	 * @return     Value of the shape function
	 */
	double quad(double x, double y, double x1, double x2, double y1, double y2, 
		unsigned int idx, valtype_e vt);

	/**
	 * @brief      Check to see if the input parameters of a 2-D triangular element is valid.
	 *
	 * @param[in]  x     independent variable x
	 * @param[in]  y     independent variable y
	 * @param[in]  x1    x coordinate of the first node of the triangle
	 * @param[in]  x2    x coordinate of the second node of the triangle
	 * @param[in]  x3    x coordinate of the third node of the triangle
	 * @param[in]  y1    y coordinate of the first node of the triangle
	 * @param[in]  y2    y coordinate of the second node of the triangle
	 * @param[in]  y3    y coordinate of the third node of the triangle
	 *
	 * @return     true if the inputs are valid.
	 */
	bool triangle_element_is_valid(double x, double y, double x1, double x2, double x3, 
		double y1, double y2, double y3);

	/**
	 * @brief      Shape function defined on a 2D triangle.
	 * 
	 * This function has three nodes on the corners, respectively.
	 * index of the nodes are as show
	 * 
	 * y
	 * |
	 * |
	 * 2
	 * | \
	 * |   \
	 * |     \
	 * |       \
	 * |         \
	 * |           \
	 * 0-------------1--->x
	 *
	 * @param[in]  x     independent variable x
	 * @param[in]  y     independent variable y
	 * @param[in]  x1    x coordinate of the first node of the triangle
	 * @param[in]  x2    x coordinate of the second node of the triangle
	 * @param[in]  x3    x coordinate of the third node of the triangle
	 * @param[in]  y1    y coordinate of the first node of the triangle
	 * @param[in]  y2    y coordinate of the second node of the triangle
	 * @param[in]  y3    y coordinate of the third node of the triangle
	 * @param[in]  idx   Node index
	 * @param[in]  vt    Calculation target. Could be Value, Dx or Dy
	 *
	 * @return     Value of the shape function
	 */
	double triangle(double x, double y, double x1, double x2, double x3, 
		double y1, double y2, double y3, unsigned int idx, valtype_e vt);

	/**
	 * @brief      Check to see if the input parameters of a 3-D cubic element is valid.
	 *
	 * @param[in]  x     independent variable x
	 * @param[in]  y     independent variable y
	 * @param[in]  z     independent variable z
	 * @param[in]  x1    low bound of x
	 * @param[in]  x2    high bound of x
	 * @param[in]  y1    low bound of y
	 * @param[in]  y2    high bound of y
	 * @param[in]  z1    low bound of z
	 * @param[in]  z2    high bound of z
	 *
	 * @return     true if the inputs are valid.
	 */
	bool cube_element_is_valid(double x, double y, double z, double x1, double x2, 
		double y1, double y2, double z1, double z2);

	/**
	 * @brief      Shape function defined on a 3D cube.
	 * 
	 * This function has three nodes on the corners, respectively.
	 * index of the nodes are as show
	 * z
	 * |
	 * |   7----------6
	 * |  /|         /|
	 * | / |        / |
	 * |/  |       /  |
	 * 4----------5   |
	 * |   3------|---2
	 * |  /       |  /
	 * | /        | /
	 * 0----------1--------> x
	 *
	 * @param[in]  x     independent variable x
	 * @param[in]  y     independent variable y
	 * @param[in]  z     independent variable z
	 * @param[in]  x1    low bound of x
	 * @param[in]  x2    high bound of x
	 * @param[in]  y1    low bound of y
	 * @param[in]  y2    high bound of y
	 * @param[in]  z1    low bound of z
	 * @param[in]  z2    high bound of z
	 * @param[in]  idx   Node index
	 * @param[in]  vt    Calculation target. Could be Value, Dx, Dy or Dz
	 *
	 * @return     Value of the shape function
	 */
	double cube(double x, double y, double z, double x1, double x2, 
		double y1, double y2, double z1, double z2, unsigned int idx, valtype_e vt);

	/**
	 * @brief      Check to see if the input parameters of a 3-D tetrahedral element is valid.
	 *
	 * @param[in]  x     independent variable x
	 * @param[in]  y     independent variable y
	 * @param[in]  z     independent variable z
	 * @param[in]  x1    x coordinate of the first node of the tetrahedron
	 * @param[in]  x2    x coordinate of the second node of the tetrahedron
	 * @param[in]  x3    x coordinate of the third node of the tetrahedron
	 * @param[in]  x4    x coordinate of the fourth node of the tetrahedron
	 * @param[in]  y1    y coordinate of the first node of the tetrahedron
	 * @param[in]  y2    y coordinate of the second node of the tetrahedron
	 * @param[in]  y3    y coordinate of the third node of the tetrahedron
	 * @param[in]  y4    y coordinate of the fourth node of the tetrahedron
	 * @param[in]  z1    z coordinate of the first node of the tetrahedron
	 * @param[in]  z2    z coordinate of the second node of the tetrahedron
	 * @param[in]  z3    z coordinate of the third node of the tetrahedron
	 * @param[in]  z4    z coordinate of the fourth node of the tetrahedron
	 *
	 * @return     true if the inputs are valid.
	 */
	//bool tetrahedron_element_is_valid(double x, double y, double z, double x1, double x2, double x3, double x4, 
	//	double y1, double y2, double y3, double y4, double z1, double z2, double z3, double z4);

	/**
	 * @brief      Shape function defined on a 3D tetrahedron.
	 * 
	 * This function has four nodes on the corners, respectively.
	 * index of the nodes are as show
	 * q
	 * |
	 * 3
	 * |\
	 * |  \
	 * |	\
	 * |	  \
	 * |		\
	 * |		  \
	 * 0------------2---> v
	 * |           /
	 * |         /
	 * |       /
	 * |     /
	 * |   /
	 * | /
	 * |/
	 * 1
	 * |
	 * u
	 *
	 * @param[in]  x     independent variable x
	 * @param[in]  y     independent variable y
	 * @param[in]  z     independent variable z
	 * @param[in]  x1    x coordinate of the first node of the tetrahedron
	 * @param[in]  x2    x coordinate of the second node of the tetrahedron
	 * @param[in]  x3    x coordinate of the third node of the tetrahedron
	 * @param[in]  x4    x coordinate of the fourth node of the tetrahedron
	 * @param[in]  y1    y coordinate of the first node of the tetrahedron
	 * @param[in]  y2    y coordinate of the second node of the tetrahedron
	 * @param[in]  y3    y coordinate of the third node of the tetrahedron
	 * @param[in]  y4    y coordinate of the fourth node of the tetrahedron
	 * @param[in]  z1    z coordinate of the first node of the tetrahedron
	 * @param[in]  z2    z coordinate of the second node of the tetrahedron
	 * @param[in]  z3    z coordinate of the third node of the tetrahedron
	 * @param[in]  z4    z coordinate of the fourth node of the tetrahedron
	 * @param[in]  idx   Node index
	 * @param[in]  vt    Calculation target. Could be Value, Dx, Dy or Dz
	 *
	 * @return     Value of the shape function
	 */
	double tetrahedron(double x, double y, double z, double x1, double x2, double x3, double x4, 
		double y1, double y2, double y3, double y4, double z1, double z2, double z3, double z4,  
		unsigned int idx, valtype_e vt);
};

}

#endif // _NFEM_SF_H